module MCV_mod

  use constants_mod
  use namelist_mod
  use boundary_mod, only: pbc_idx

  implicit none
  
  private
  public :: compute_flux_jacobian_MCV
  public :: ode_rhs_mcv

  real(r8), parameter :: coef_L(1:4) = [ 1.0_r8, -6.0_r8, 3.0_r8,  2.0_r8] 
  real(r8), parameter :: coef_R(1:4) = [-2.0_r8, -3.0_r8, 6.0_r8, -1.0_r8]
  
  contains

  subroutine compute_flux_jacobian_MCV(q, dt, jac)
    
    real(r8), intent(in) :: q(1-nhalo:nx+nhalo)
    real(r8), intent(in) :: dt
    real(r8), intent(out) :: jac(1:nx, 1:nx)
    integer :: i, idx_m3, idx_m2, idx_m1, idx_p1, idx_p2, idx_p3
    real(r8) :: alpha

    alpha = abs(u)
    jac = 0.0_r8

    ! Odd points (face-centered)
    do i = 1, nx, 2
      ! Get periodic indices according to exact PBC rules
      idx_m2 = pbc_idx(i-2, nx)
      idx_m1 = pbc_idx(i-1, nx)
      idx_p1 = pbc_idx(i+1, nx)
      idx_p2 = pbc_idx(i+2, nx)

      ! Jacobian entries for odd points
      jac(i, idx_m2) = -0.5_r8 * (          u +          alpha) / (3.0_r8 * dx) ! d(f_xi)/dq_{i-2}
      jac(i, idx_m1) = -0.5_r8 * (-8.0_r8 * u - 4.0_r8 * alpha) / (3.0_r8 * dx) ! d(f_xi)/dq_{i-1}
      jac(i, i     ) = -0.5_r8 * (              6.0_r8 * alpha) / (3.0_r8 * dx) ! d(f_xi)/dq_i
      jac(i, idx_p1) = -0.5_r8 * ( 8.0_r8 * u - 4.0_r8 * alpha) / (3.0_r8 * dx) ! d(f_xi)/dq_{i+1}
      jac(i, idx_p2) = -0.5_r8 * (         -u +          alpha) / (3.0_r8 * dx) ! d(f_xi)/dq_{i+2}
    end do

    ! Even points (center-centered)
    do i = 2, nx-1, 2
      ! Get periodic indices according to exact PBC rules
      idx_m3 = pbc_idx(i-3, nx)
      idx_m2 = pbc_idx(i-2, nx)
      idx_m1 = pbc_idx(i-1, nx)
      idx_p1 = pbc_idx(i+1, nx)
      idx_p2 = pbc_idx(i+2, nx)
      idx_p3 = pbc_idx(i+3, nx)
      ! j = i -3
      jac(i, idx_m3) =                   (          u +          alpha) / (24.0_r8 * dx)
      jac(i, idx_m2) =                   (-8.0_r8 * u - 4.0_r8 * alpha) / (24.0_r8 * dx)
      jac(i, idx_m1) = 1.5_r8 * u / dx + (          u + 7.0_r8 * alpha) / (24.0_r8 * dx)
      jac(i, i     ) =                                  -8.0_r8 * alpha / (24.0_r8 * dx)
      jac(i, idx_p1) = -1.5_r8 * u / dx + (        -u + 7.0_r8 * alpha) / (24.0_r8 * dx)
      jac(i, idx_p2) =                    (8.0_r8 * u - 4.0_r8 * alpha) / (24.0_r8 * dx)
      jac(i, idx_p3) =                    (        -u +          alpha) / (24.0_r8 * dx)
    end do

  end subroutine compute_flux_jacobian_MCV

  function ode_rhs_mcv(var, dt) result(rhs)

    real(r8), intent(in) :: var(1-nhalo:nx+nhalo)
    real(r8), intent(in) :: dt
    real(r8) :: rhs(1:nx)
    integer :: i
    real(r8) :: alpha

    alpha = abs(u)

    do i = 1, nx, 2
      rhs(i) = -0.5_r8 * (u * (dot_product(coef_L, var(i-2:i+1)) + dot_product(coef_R, var(i-1:i+2))) - &
                      alpha * (dot_product(coef_R, var(i-1:i+2)) - dot_product(coef_L, var(i-2:i+1)))) / (3.0_r8 * dx)
    end do
    
    do i = 2, nx-1, 2
      rhs(i) = -1.5_r8 * u * (var(i+1) - var(i-1)) / dx + 0.25_r8 * (-rhs(i+1) - rhs(i-1))
    end do

  end function ode_rhs_mcv

end module MCV_mod